function f = bifsstheta(p0, gamma0, thetaLower, thetaUpper)

g = 9.8;
m = 100;
L = 2*pi;
h = (thetaUpper-thetaLower)/m;
if (gamma0 < 0)
	for i = 1:m
		theta(i) = (i-1)*h+ thetaLower;
		f(i) = (4*pi^2*sqrt(-2*gamma0)/(3*L^2))*[(p0*theta(i)/100 - p0)^3] + ...
			(-2*gamma0)^(1.5)*(p0*theta(i)/100 - p0) - g*p0^2;
	end
else 
	for i = 1:m
		theta(i) = (i-1)*h+ thetaLower;
		f(i) = (4*pi^2*sqrt(2*gamma0)/(3*L^2))*[(-p0)^3-...
			(p0*theta(i)/100 - p0)^3] + ...
			(2*gamma0)^(1.5)*(-p0*theta(i)/100) - g*p0^2;
	end
end
plot(theta, f);
xlabel('\theta', 'FontSize', 14);
title({['p0 = ', num2str(p0), ', \gamma_0 = ', gamma0]});
end


